%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Explanation of what the file does
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% This function solves for the entire era (long timespan of T years)
% equilibrium of dots of price, actual labors, and wages, and levels of
% revenues and trade shares given labor supply and initial conditions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Inputs needed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% 1. Deficit matrix (DefEra), size I*T
% 2. Initial income by sector (labincsec), size I*S
% 3. Initial revenues (revbysec), size I*S
% 4. Initial trade shares (tradesharedat), size I*I*S
% 5. Dot trade costs (TAUmatC), size I*I*S*T
% 6. Dot technologies (TECmatC), size I*S*T
% 7. Labor shares (labshares), size I*S
% 8. IO matrix (inoutmat), size I*S*S
% 9. Trade elasticities (Smat), size 1*S
% 10. Final good consumption shares (Amat), size I*S
% 11. Labor supply matrix (EllEra), size I*(S+1)*T
% 12. Delta DNWR matrix (DeltaMat), size I*S
% 13. Parameter for updating AL algorithm (lambda), scalar
% 14. Tolerance parameter (tol), scalar
% 15. Agglomeration w.r.t. regional employment (AECR), flexible row vector
% 16. Agglome w.r.t. region-sector employment (AECRS), flexible row vector
% 17. Agglomeration parameter (aeatw), scalar, 1-> everyone, 2-> just USA
% 18. Place where non-USA regions start (M), scalar
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Outputs produced
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% 1. Wages dot (Wmat), size I*S*(T-1)
% 2. Labor dot (Lmat), size I*S*(T-1)
% 3. Prices dot (PHma), size I*S*(T-1)
% 4. New revenue levels (RELE), size I*S*T
% 5. New trade share levels (TSLE), size I*I*S*T
% 6. Income by sector (Ytm1), size I*S*T
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Matlab functions invoked
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% 1. NEBA_PeriodTatoSolver

function [Wmat,Lmat,PHma,RELE,TSLE,Ytm1]= NEBA_Block2(DefEra,labincsec,...
revbysec,tradesharedat,TAUmatC,TECmatC,labshares,inoutmat,Smat,...
Amat,EllEra,DeltaMat,lambda,tol,AECR,AECRS,aeatw,M)

% Define some useful variables

EllRel              = EllEra(:,2:end,:);
I                   = size(EllRel,1);
S                   = size(EllRel,2);
T                   = size(EllRel,3);

% Declare matrices of results

Wmat                = zeros(I,S,T-1);
Lmat                = zeros(I,S,T-1);
PHma                = zeros(I,S,T-1);

% Declare matrices of past levels

Ytm1                = zeros(I,S,T);
RELE                = zeros(I,S,T);
TSLE                = zeros(I,I,S,T);

% Fill some matrices in levels for the initial period (from the data)

Ytm1(:,:,1)         = labincsec;
RELE(:,:,1)         = revbysec;
TSLE(:,:,:,1)       = tradesharedat;

% Now run the time loop

for t=1:T-1

    Deficits        = DefEra(:,t+1);
    PastYnom        = Ytm1(:,:,t);
    PastSha         = TSLE(:,:,:,t);
    
    if t==1
        PastCumLhat      = ones(I,S);
        TECmatrix        = TECmatC(:,:,t);
    else
        PastCumLhat      = prod(Lmat(:,:,1:t-1),3);
        CumLhat(:,:,1)   = ones(I,S);
        CumLhat(:,:,2:t) = cumprod(Lmat(:,:,1:t-1),3);
        LLev             = CumLhat.*repmat(EllRel(:,:,1),1,1,t);
        LLevRegion       = sum(LLev,2);
        LRegDot          = LLevRegion(:,:,2:end)./LLevRegion(:,:,1:end-1);
        LRegDotrep       = repmat(LRegDot,1,S,1); 
        sizem            = min(t-1,size(AECR,2));
        Relexponent      = repmat(permute(AECR( end-sizem+1:end),[1 3 2]),I,S,1);
        RaLRegDotrep     = LRegDotrep(:,:,end-sizem+1:end).^Relexponent;
        RelLmat          = Lmat(:,:,1:t-1);
        Relexponent2     = repmat(permute(AECRS(end-sizem+1:end),[1 3 2]),I,S,1);
        RaLDot           = RelLmat(:,:,end-sizem+1:end).^Relexponent2;
        TECmatrix        = TECmatC(:,:,t) .* prod(RaLRegDotrep,3) .* prod(RaLDot,3);
    end
    
    if aeatw==2
        TECmatrix(M+1:end,:) = TECmatC(M+1:end,:,t);
    end
    
    TAUmatrix       = TAUmatC(:,:,:,t);

    % Find the equilibrium for period t using a specially designed function

    Llim            = (EllRel(:,:,t+1)./EllRel(:,:,1))./PastCumLhat;
    maxiter         = 15000;
    [EqW,EqL,EqP,EqS,EqR,~] = NEBA_PeriodTatoSolver(PastSha,TAUmatrix,...
    TECmatrix,labshares,inoutmat,Smat,Amat,PastYnom,Llim,DeltaMat,...
    Deficits,maxiter,tol,lambda);

    % Update some matrices like wages, employment, prices, revenues, etc

    Wmat(:,:,t)     = EqW;
    Lmat(:,:,t)     = EqL;
    PHma(:,:,t)     = EqP;
    Ytm1(:,:,t+1)   = EqW.*EqL.*PastYnom;
    RELE(:,:,t+1)   = EqR;
    TSLE(:,:,:,t+1) = EqS;
        
end

end
